In many practical applications it is required to obtain and analyze three-dimensional (3D) arrays of data representing one or more physical properties at regular grid positions within the interior of solid bodies. Typically, the data is obtained by methods such as computer axial tomography (CAT), x-ray scanning, nuclear magnetic resonance (NMR) imaging systems, or by other non-intrusive mechanisms such as ultrasound scanning, positron emission tomography (PET), emission computer tomography (ECT) and multimodality imaging (MMI). Each of these techniques produces a planar, grid-like array of values for a continuous sequence of planes or slices of the solid object, thus providing a three-dimensional array of values. A typical application is in medicine where the solid object can be a human body or a part of it and the objective is to detect and model certain internal structures for diagnostic and surgical planning purposes. The methods to which the present invention relates are, however, equally applicable to other natural or artificial bodies and fields of study.
Examples of the physical properties measured at regular 3D positions include the coefficient of x-ray absorption in the CAT case or, for NMR imaging, the spin-spin or the spin-lattice relaxation parameters. In each case, the measured physical values reflect variations in the composition, density or structural characteristics of the underlying physical objects, providing knowledge about their internal structure. The 3D data array typically corresponds to a plurality of measurements of the given physical property, taken at regular intervals along an appropriately selected system of 3D (X, Y, Z) physical coordinates. Each measurement V.sub.XYZ of the physical property is thus associated with a specific array index position in a cubic lattice within the volume in study. A cubically adjacent set of eight such positions defines an unit cubic volume of the body under study and is known in the art as a voxel.
It is further known in the prior art to utilize such three-dimensional arrays of interior physical values to generate visual images of the structures within the body on a display system. In order to display two-dimensional (2D) images of such 3D interior structures, however, it is necessary to locate the position of the surface of such structure within the array of physical values. This is accomplished by comparing the array values to a threshold value or to a range of threshold values which correspond to the physical property values normally associated with that surface.
For example, in the case of the human body, bones, brain or any other tissues can be characterized by the known range of their density values which is compared to the measured and stored 3D array values in order to separate out a surface model of their structure. The displayed visual images can next be used to assist the medical practitioner in providing a diagnosis or in the planning of complex surgical procedures.
Once the location of the surface of an internal structure is determined, this surface is typically illuminated and shaded so as to give the human eye the impression of a natural shape and spatial disposition of the surface when it is displayed. Normally, to provide such shading the angular direction of a vector orthogonal to the surface at each point on the surface (the surface normal) is compared to the viewing angle of the observer. The intensity of shading can be adjusted to be proportional to the difference between these angles.
There are several problems associated with the surface rendering of structures in the interior of a solid body. First, it is necessary to accurately determine which of the voxel elements of the 3D array have been crossed by the surface of interest and at which physical location. Next, it is important to provide a mechanism for displaying the internal structure of a three-dimensional body on the two-dimensional computer or TV screen. Another important problem is that data acquired from physical data collections systems is often noisy and is not capable of being directly displayed prior to some kind of a filtering operation.
Several solutions have been proposed in the past to address these problems. An approach that is frequently used is the one described in U.S. Pat. No. 4,710,876 to Cline et al. on the so-called "marching cubes". In this method, the surface segment intersecting a voxel is approximated by one of a limited number of standardized plane polygonal surfaces intersecting the voxel. One particular standardized surface tiling is selected from a library of available surfaces by a vector which represents the binary differences between the threshold value for the given surface of interest and the eight voxel vertex values. The surface-to-voxel intersection coordinates, as well as the vector normal to the surface for each such standardized polygonal surface tiling can be obtained by table look-up techniques. The final surface is assembled as a mosaic using all the standardized triangular tiles.
The marching cubes method has found significant acceptance in the industry because it is relatively simple to implement and is fast. Several modifications have been proposed since 1987. For instance, in U.S. Pat. No. 4,729,098 granted in March 1988 to Cline et al. a variation of the marching cubes methods is proposed, in which non-linear interpolation is used to more accurately locate the coordinates of the standardized polygons.
Another known method for approximating the surface of an internal structure is the so-called "dividing cubes" method, disclosed in U.S. Pat. No. 4,719,585 to Cline et al. In this method, the values at the vertices of the voxel are used to interpolate in three dimensions values at regularly positioned intra-voxel sub-grid locations. The interpolated sub-grid values can then be used to locate the surface position and calculate the normal vector more accurately.
In order to distinguish between different internal structures having the same or similar physical property and thus having close or similar values on the three-dimensional grid, U.S. Pat. No. 4,751,643 to Lorensen discloses a technique for labeling surfaces with similar property values and using adjacency criteria in respect to a "seed location" in the particular structure of interest to segregate the desired surface from all of the labeled surfaces.
U.S. Pat. No. 4,791,567 to Cline et al. discloses another technique of segregating similar internal structures by determining connectivity from adjacency information. U.S. Pat. No. 4,821,213 to Cline et al. discloses yet another technique for differentiating internal structures in which a linear path is made through the data array to locate and label all different structures along the scan line by counting structure interfaces.
In the methods described above, voxels within the solid object under study are typically traversed sequentially and progressively, so that every voxel within the data space is checked for a crossing by the surface of interest, and is subsequently assigned to belong either to the outside or to the inside of the surface. Because of the large volumes of data to be processed, the methods are computationally intensive and relatively slow. An improvement in this respect is described in U.S. Pat. No. 5,166,876 to Cline, issued in November 1992. The disclosed method significantly speeds up the rendering of a surface within a boundary by only processing voxels on the surface. This approach saves significant computation time but still uses the same case table and the limited set of approximating surfaces which were described in the original marching cubes patent.
As noted above, the marching cubes surface rendering method is widely accepted but has certain deficiencies. These deficiencies are generally related to the fact that the fifteen canonical cases used in the surface segments' approximation do not cover the entire spectrum of possible tiling combinations. This results in occasional approximation errors, occurring at lower levels of resolution, specifically if a voxel face is intersected by the surface twice resulting in four crossings ("hits") with the voxel edge. The errors manifest themselves in incorrect surface rendering, visually appearing as unexpected spots, tears and missing tile elements. The problems are aggravated if the collection data is noisy and/or with low contrast.
A number of corrections to the marching cube algorithm have been offered in the past. Notable is a paper by Neilson and Hamann, "The Asymptotic Decider: Resolving the Ambiguity in Marching Cubes," in the Proceedings of Visualization 1991, pp. 83-91. This and other publications show that the case table approach produces inconsistent results in some cases. Thus, no topologically correct algorithm for isosurface tiling has been offered so far. Accordingly, reliable measurements of internal structure parameters, such as surface area and volume, are also not available.
The display of 3D graphic images has been principally driven by the goals and directions of computer aided design (CAD) and computer aided manufacturing (CAM) systems. Systems accepting data in several basic formats have been developed for displaying solid bodies and manipulating images in various ways to create models of parts being manufactured, or internal organs being monitored.
The present systems for surface rendering and display are, however, still incapable of meeting the requirements for a fast and truly interactive display of continuously monitored 3D data. The complexity of the task can be better appreciated by noting that a typical surface rendering operation involves processing on the order of 512.times.512.times.512 pixels. In the presence of noise the task becomes even more complicated because of the necessity for additional noise filtering, smoothing and data interpolation, which takes processing time and often results in the loss of important structural details.
It is thus desirable to provide a method for rendering topologically correct surface approximations of arbitrary internal structures within a volume of interest and a system for the display and analysis of such surface structures. It is also desirable to provide high quality, interactive display of the surface structures of interest from a set of collected data and minimize the effects of possible noise interference. Additionally, it is desirable to provide an efficient implementation of the method which facilitates the use of specialized hardware for real time data processing.